pSERG CIRCADIAN DISTRIBUTION OF STATUS EPILEPTICUS
This analysis evaluates the circadian distribution of refractory status epilepticus episodes
Install packages needed for analysis
# install.packages("plyr")
library(plyr)
# install.packages("gdata")
library(gdata)
## gdata: Unable to locate valid perl interpreter
## gdata:
## gdata: read.xls() will be unable to read Excel XLS and XLSX files
## gdata: unless the 'perl=' argument is used to specify the location
## gdata: of a valid perl intrpreter.
## gdata:
## gdata: (To avoid display of this message in the future, please
## gdata: ensure perl is installed and available on the executable
## gdata: search path.)
## gdata: Unable to load perl libaries needed by read.xls()
## gdata: to support 'XLX' (Excel 97-2004) files.
##
## gdata: Unable to load perl libaries needed by read.xls()
## gdata: to support 'XLSX' (Excel 2007+) files.
##
## gdata: Run the function 'installXLSXsupport()'
## gdata: to automatically download and install the perl
## gdata: libaries needed to support Excel XLS and XLSX formats.
##
## Attaching package: 'gdata'
## The following object is masked from 'package:stats':
##
## nobs
## The following object is masked from 'package:utils':
##
## object.size
## The following object is masked from 'package:base':
##
## startsWith
# install.packages("gmodels")
library(gmodels)
# install.packages("cosinor")
library(cosinor)
## Warning: package 'cosinor' was built under R version 3.5.3
# install.packages("ggplot2")
library(ggplot2)
# install.packages("plotly")
library(plotly)
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following objects are masked from 'package:plyr':
##
## arrange, mutate, rename, summarise
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
Load the database
## Load the data
pSERGcircadian <- read.csv("F:\\PROJECTS\\pSERG Circadian\\Data\\pSERG.csv")
Data cleaning and preparation of variables for analysis
## Keep patients only if they are a refractory SE case (not a control)
pSERGcircadian <- pSERGcircadian[which(pSERGcircadian$SE_GROUP == 'refractory_case'), ]
pSERGcircadian$SE_GROUP <- droplevels(pSERGcircadian$SE_GROUP)
## Transform date of status epilepticus into date format
pSERGcircadian$DATESEIZURE <- as.Date(pSERGcircadian$DATESEIZURE, format = "%m/%d/%Y")
## Order by date of status epilepticus
pSERGcircadian <- pSERGcircadian[order(pSERGcircadian$PATIENT_LABEL, pSERGcircadian$DATESEIZURE), ]
## Delete duplicate episodes from the same patient
pSERGcircadian <- pSERGcircadian[!duplicated(pSERGcircadian$PATIENT_LABEL), ]
#################Create variable hour of SE######################
## Transform time of hour to numeric and delete cases with unknown hour of SE
pSERGcircadian$TIMESEIZURE_HOURS <- as.numeric(as.character(pSERGcircadian$TIMESEIZURE_HOURS))
## Warning: NAs introduced by coercion
pSERGcircadian$TIMESEIZURE_HOURS <- ifelse(pSERGcircadian$TIMESEIZURE_HOURS>99, pSERGcircadian$TIMESEIZURE_HOURS%/%100, pSERGcircadian$TIMESEIZURE_HOURS)
## Delete patients with no numeric values for hour of seizure onset or no information
pSERGcircadian <- pSERGcircadian[!is.na(pSERGcircadian$TIMESEIZURE_HOURS), ]
## Rename variable
pSERGcircadian$hourofSE <- pSERGcircadian$TIMESEIZURE_HOURS
## Rename 0 to 24
pSERGcircadian$hourofSE[pSERGcircadian$hourofSE == 0] <- 24
## Exact hour
pSERGcircadian$TIMESEIZURE_MINUTES <- as.numeric(as.character(pSERGcircadian$TIMESEIZURE_MINUTES))
## Warning: NAs introduced by coercion
pSERGcircadian$time <- pSERGcircadian$TIMESEIZURE_HOURS + (pSERGcircadian$TIMESEIZURE_MINUTES / 60)
#######################Transform age into numeric and delete patients with unknown age####################
# Transform age into numeric
pSERGcircadian$G_T_STTS_PLPTCS_EPISODE_MONTHS <- as.numeric(as.character(pSERGcircadian$G_T_STTS_PLPTCS_EPISODE_MONTHS))
## Warning: NAs introduced by coercion
pSERGcircadian$G_T_STTS_PLPTCUS_EPISODE_YEARS <- as.numeric(as.character(pSERGcircadian$G_T_STTS_PLPTCUS_EPISODE_YEARS))
## Warning: NAs introduced by coercion
pSERGcircadian$ageyears <- pSERGcircadian$G_T_STTS_PLPTCUS_EPISODE_YEARS + (pSERGcircadian$G_T_STTS_PLPTCS_EPISODE_MONTHS/12)
pSERGcircadian <- pSERGcircadian[complete.cases(pSERGcircadian[ , "ageyears"]), ]
##########################Delete patients with unknown sex################################################
pSERGcircadian <- pSERGcircadian[which(pSERGcircadian$SEX == "male" | pSERGcircadian$SEX == "female"), ]
#######################Transform variables for analysis###############
## Create variable delay
pSERGcircadian$delay[grepl("delay", pSERGcircadian$PAST)] <- 1
pSERGcircadian$delay[!grepl("delay", pSERGcircadian$PAST)] <- 0
## Create variable palsy
pSERGcircadian$palsy[grepl("palsy", pSERGcircadian$PAST)] <- 1
pSERGcircadian$palsy[!grepl("palsy", pSERGcircadian$PAST)] <- 0
## Create variable febrile
pSERGcircadian$febrile[grepl("febrile", pSERGcircadian$PAST)] <- 1
pSERGcircadian$febrile[!grepl("febrile", pSERGcircadian$PAST)] <- 0
## Create variable prior epilepsy
pSERGcircadian$priorepilepsy[grepl("epi", pSERGcircadian$PAST)] <- 1
pSERGcircadian$priorepilepsy[!grepl("epi", pSERGcircadian$PAST)] <- 0
## Create variable status
pSERGcircadian$status[grepl("status", pSERGcircadian$PAST)] <- 1
pSERGcircadian$status[!grepl("status", pSERGcircadian$PAST)] <- 0
## Create variable none
pSERGcircadian$none[grepl("none", pSERGcircadian$PAST)] <- 1
pSERGcircadian$none[!grepl("none", pSERGcircadian$PAST)] <- 0
Demographics
############################DEMOGRAPHICS#######################
## Gender
CrossTable(pSERGcircadian$SEX)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 368
##
##
## | female | male |
## |-----------|-----------|
## | 155 | 213 |
## | 0.421 | 0.579 |
## |-----------|-----------|
##
##
##
##
## Age
nobs(pSERGcircadian$ageyears)
## [1] 368
summary(pSERGcircadian$ageyears)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.08333 1.32146 4.16667 5.94850 9.68540 20.74167
sd(pSERGcircadian$ageyears, na.rm = TRUE)
## [1] 5.27164
# Patients 18 year-old or older
pSERGcircadian$eighteenormore <- as.factor(ifelse(pSERGcircadian$ageyears >= 18, 1, 0))
CrossTable(pSERGcircadian$eighteenormore)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 368
##
##
## | 0 | 1 |
## |-----------|-----------|
## | 361 | 7 |
## | 0.981 | 0.019 |
## |-----------|-----------|
##
##
##
##
## Race
CrossTable(pSERGcircadian$RACE)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 368
##
##
## | american_indian_alaska_native | arabic | asian | black_or_african_american | native_hawaiian_or_pacific_islander |
## |-------------------------------------|-------------------------------------|-------------------------------------|-------------------------------------|-------------------------------------|
## | 1 | 10 | 11 | 71 | 1 |
## | 0.003 | 0.027 | 0.030 | 0.193 | 0.003 |
## |-------------------------------------|-------------------------------------|-------------------------------------|-------------------------------------|-------------------------------------|
##
##
## | not_reported | unknown | white |
## |-------------------------------------|-------------------------------------|-------------------------------------|
## | 15 | 22 | 237 |
## | 0.041 | 0.060 | 0.644 |
## |-------------------------------------|-------------------------------------|-------------------------------------|
##
##
##
##
## Ethnicity
CrossTable(pSERGcircadian$ETHNICITY)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 368
##
##
## | hispanic_or_latino | not_hispanic_or_latino | not_reported | unknown |
## |------------------------|------------------------|------------------------|------------------------|
## | 66 | 269 | 21 | 12 |
## | 0.179 | 0.731 | 0.057 | 0.033 |
## |------------------------|------------------------|------------------------|------------------------|
##
##
##
##
## Delay
CrossTable(pSERGcircadian$delay)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 368
##
##
## | 0 | 1 |
## |-----------|-----------|
## | 182 | 186 |
## | 0.495 | 0.505 |
## |-----------|-----------|
##
##
##
##
## Palsy
CrossTable(pSERGcircadian$palsy)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 368
##
##
## | 0 | 1 |
## |-----------|-----------|
## | 329 | 39 |
## | 0.894 | 0.106 |
## |-----------|-----------|
##
##
##
##
## Febrile
CrossTable(pSERGcircadian$febrile)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 368
##
##
## | 0 | 1 |
## |-----------|-----------|
## | 328 | 40 |
## | 0.891 | 0.109 |
## |-----------|-----------|
##
##
##
##
## Prior epilepsy
CrossTable(pSERGcircadian$priorepilepsy)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 368
##
##
## | 0 | 1 |
## |-----------|-----------|
## | 186 | 182 |
## | 0.505 | 0.495 |
## |-----------|-----------|
##
##
##
##
## Status
CrossTable(pSERGcircadian$status)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 368
##
##
## | 0 | 1 |
## |-----------|-----------|
## | 299 | 69 |
## | 0.812 | 0.188 |
## |-----------|-----------|
##
##
##
##
## None
CrossTable(pSERGcircadian$none)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 368
##
##
## | 0 | 1 |
## |-----------|-----------|
## | 248 | 120 |
## | 0.674 | 0.326 |
## |-----------|-----------|
##
##
##
##
## Transform CONVULSIVEDURATION into numeric
pSERGcircadian$CONVULSIVEDURATION <- as.numeric(as.character(pSERGcircadian$CONVULSIVEDURATION))
## Warning: NAs introduced by coercion
pSERGcircadian$CONVULSIVEmin <- pSERGcircadian$CONVULSIVEDURATION
pSERGcircadian$CONVULSIVEhr <- pSERGcircadian$CONVULSIVEDURATION*60
pSERGcircadian$convulsivedurationmin <- ifelse(pSERGcircadian$CONVULSIVEDURATIONUNITS == "min", pSERGcircadian$CONVULSIVEmin, pSERGcircadian$CONVULSIVEhr)
## Duration of SE
nobs(pSERGcircadian$convulsivedurationmin)
## [1] 355
summary(pSERGcircadian$convulsivedurationmin)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0 60 122 1925 283 172800 13
sd(pSERGcircadian$convulsivedurationmin, na.rm = TRUE)
## [1] 12308.02
## Hospital onset
CrossTable(pSERGcircadian$HOSPITALONSET)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 368
##
##
## | no | yes |
## |-----------|-----------|
## | 258 | 110 |
## | 0.701 | 0.299 |
## |-----------|-----------|
##
##
##
##
## Type of status epilepticus
CrossTable(pSERGcircadian$TYPESTATUS)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 368
##
##
## | continuous | intermittent |
## |--------------|--------------|
## | 116 | 252 |
## | 0.315 | 0.685 |
## |--------------|--------------|
##
##
##
##
Figures of the distribution of status epilepticus episodes
## Create dummy variable for events of SE
pSERGcircadian$events <- 1
################Every hour###############
## Draw the number of SE episodes in a plot with stacked units
distributionhour <- ggplot(pSERGcircadian, aes(x = hourofSE, y = events)) + geom_bar(stat = 'identity', color = "#0076c0", fill = "#67771a") +
scale_x_discrete(limits = c(1,24)) +
theme(panel.background = element_rect(fill = "white"),
axis.text = element_text(size = 18, color = "black", face = "bold"),
axis.title = element_text(size = 24, color = "black", face = "bold"),
axis.ticks.length = unit(0.4, "cm"), axis.ticks = element_line(size = 1.2)) +
scale_x_continuous(breaks = c(4, 8, 12, 16, 20, 24)) +
scale_y_continuous(breaks = c(5, 10, 15, 20, 25)) +
labs(x= "Time of day (hour)", y= "Number of rSE episodes")
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.
distributionhour

# Interactive graph
ggplotly(distributionhour)
################Every 2 hours###############
# Create variable every 2 hours
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 1] <- "1-2"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 2] <- "1-2"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 3] <- "3-4"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 4] <- "3-4"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 5] <- "5-6"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 6] <- "5-6"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 7] <- "7-8"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 8] <- "7-8"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 9] <- "9-10"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 10] <- "9-10"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 11] <- "11-12"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 12] <- "11-12"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 13] <- "13-14"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 14] <- "13-14"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 15] <- "15-16"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 16] <- "15-16"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 17] <- "17-18"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 18] <- "17-18"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 19] <- "19-20"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 20] <- "19-20"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 21] <- "21-22"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 22] <- "21-22"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 23] <- "23-24"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 24] <- "23-24"
# Give this variable an order
pSERGcircadian$groups2 <- factor(pSERGcircadian$groups2, c("1-2", "3-4","5-6", "7-8","9-10", "11-12","13-14", "15-16","17-18", "19-20","21-22", "23-24"))
## Draw the number of SE episodes in a plot with stacked units
distribution2hours <- ggplot(pSERGcircadian[!is.na(pSERGcircadian$groups2), ], aes(x = groups2, y = events))+geom_bar(stat ='identity', color = "#0076c0", fill = "#67771a") +
theme(panel.background = element_rect(fill = "white"),
axis.text = element_text(size = 18, color = "black", face = "bold"),
axis.title = element_text(size = 24, color = "black", face = "bold"),
axis.ticks.length = unit(0.4, "cm"), axis.ticks = element_line(size = 1.2)) +
scale_x_discrete() + scale_y_continuous(breaks = c(10, 20, 30, 40)) +
labs(x= "Time of day (hour)", y= "Number of rSE episodes")
distribution2hours

# Interactive graph
ggplotly(distribution2hours)
################Every 3 hours###############
# Create variable every 3 hours
pSERGcircadian$groups3[pSERGcircadian$hourofSE == 1] <- "1-3"
pSERGcircadian$groups3[pSERGcircadian$hourofSE == 2] <- "1-3"
pSERGcircadian$groups3[pSERGcircadian$hourofSE == 3] <- "1-3"
pSERGcircadian$groups3[pSERGcircadian$hourofSE == 4] <- "4-6"
pSERGcircadian$groups3[pSERGcircadian$hourofSE == 5] <- "4-6"
pSERGcircadian$groups3[pSERGcircadian$hourofSE == 6] <- "4-6"
pSERGcircadian$groups3[pSERGcircadian$hourofSE == 7] <- "7-9"
pSERGcircadian$groups3[pSERGcircadian$hourofSE == 8] <- "7-9"
pSERGcircadian$groups3[pSERGcircadian$hourofSE == 9] <- "7-9"
pSERGcircadian$groups3[pSERGcircadian$hourofSE == 10] <- "10-12"
pSERGcircadian$groups3[pSERGcircadian$hourofSE == 11] <- "10-12"
pSERGcircadian$groups3[pSERGcircadian$hourofSE == 12] <- "10-12"
pSERGcircadian$groups3[pSERGcircadian$hourofSE == 13] <- "13-15"
pSERGcircadian$groups3[pSERGcircadian$hourofSE == 14] <- "13-15"
pSERGcircadian$groups3[pSERGcircadian$hourofSE == 15] <- "13-15"
pSERGcircadian$groups3[pSERGcircadian$hourofSE== 16] <- "16-18"
pSERGcircadian$groups3[pSERGcircadian$hourofSE== 17] <- "16-18"
pSERGcircadian$groups3[pSERGcircadian$hourofSE== 18] <- "16-18"
pSERGcircadian$groups3[pSERGcircadian$hourofSE== 19] <- "19-21"
pSERGcircadian$groups3[pSERGcircadian$hourofSE== 20] <- "19-21"
pSERGcircadian$groups3[pSERGcircadian$hourofSE== 21] <- "19-21"
pSERGcircadian$groups3[pSERGcircadian$hourofSE== 22] <- "22-24"
pSERGcircadian$groups3[pSERGcircadian$hourofSE== 23] <- "22-24"
pSERGcircadian$groups3[pSERGcircadian$hourofSE== 24] <- "22-24"
# Give this variable an order
pSERGcircadian$groups3 <- factor(pSERGcircadian$groups3, c("1-3", "4-6", "7-9", "10-12", "13-15", "16-18", "19-21", "22-24"))
## Draw the number of SE in a plot
distribution3hours <- ggplot(pSERGcircadian[!is.na(pSERGcircadian$groups3), ], aes(x = groups3, y = events)) + geom_bar(stat = 'identity', color = "#0076c0", fill = "#67771a") +
theme(panel.background=element_rect(fill = "white"),
axis.text = element_text(size = 18, color = "black", face = "bold"),
axis.title = element_text(size = 24, color = "black", face = "bold"),
axis.ticks.length = unit(0.4, "cm"), axis.ticks = element_line(size = 1.2)) +
scale_x_discrete() + scale_y_continuous(breaks = c(10, 20, 30, 40, 50)) +
labs(x = "Time of day (hour)", y = "Number of rSE episodes")
distribution3hours

# Interactive graph
ggplotly(distribution3hours)
################Every 4 hours###############
#Create variable every 4 hours
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 1] <- "1-4"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 2] <- "1-4"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 3] <- "1-4"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 4] <- "1-4"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 5] <- "5-8"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 6] <- "5-8"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 7] <- "5-8"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 8] <- "5-8"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 9] <- "9-12"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 10] <- "9-12"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 11] <- "9-12"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 12] <- "9-12"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 13] <- "13-16"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 14] <- "13-16"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 15] <- "13-16"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 16] <- "13-16"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 17] <- "17-20"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 18] <- "17-20"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 19] <- "17-20"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 20] <- "17-20"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 21] <- "21-24"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 22] <- "21-24"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 23] <- "21-24"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 24] <- "21-24"
#Give this variable an order
pSERGcircadian$groupsfour <- factor(pSERGcircadian$groupsfour, c("1-4","5-8","9-12","13-16","17-20","21-24"))
##Draw the number of SE in a plot
distribution4hours <- ggplot(pSERGcircadian[!is.na(pSERGcircadian$groupsfour), ], aes(x = groupsfour, y = events)) + geom_bar(stat = 'identity', color = "#0076c0", fill = "#67771a") +
theme(panel.background = element_rect(fill = "white"),
axis.text = element_text(size = 18, color = "black", face = "bold"),
axis.title = element_text(size = 24, color = "black", face = "bold"),
axis.ticks.length = unit(0.4, "cm"), axis.ticks = element_line(size = 1.2)) +
scale_x_discrete() + scale_y_continuous(breaks = c(10, 20, 30, 40, 50, 60, 70)) +
labs(x = "Time of day (hour)", y = "Number of rSE episodes")
distribution4hours

# Interactive graph
ggplotly(distribution4hours)
Primary outcome
# Test (Kolmogorov-Smirnov) if the distribution of SE during the day is consistent with a uniform distribution
ks.test(pSERGcircadian$hourofSE, "punif", 0, 24)
## Warning in ks.test(pSERGcircadian$hourofSE, "punif", 0, 24): ties should
## not be present for the Kolmogorov-Smirnov test
##
## One-sample Kolmogorov-Smirnov test
##
## data: pSERGcircadian$hourofSE
## D = 0.088768, p-value = 0.006058
## alternative hypothesis: two-sided
####################Distribution of SE episodes through the day###########################
# Distribution of SE episodes if they were evenly distributed
nobs(pSERGcircadian$events)
## [1] 368
nobs(pSERGcircadian$events)/24
## [1] 15.33333
#####################################################################################
# Create another database with the circadian data
occurrences <- count(pSERGcircadian, "hourofSE")
# Make it a dataframe
circadiandata <- data.frame(occurrences)
circadiandata$numberSE <- circadiandata$freq
circadiandata <- circadiandata[ , c("hourofSE", "numberSE")]
#####################################################################################
# Cosinor analysis
SEfit <- cosinor.lm(numberSE ~ time(hourofSE), data = circadiandata, period = 24)
summary(SEfit)
## Raw model coefficients:
## estimate standard.error lower.CI upper.CI p.value
## (Intercept) 15.3333 0.7472 13.8688 16.7978 0.0000
## rrr -3.1912 1.0567 -5.2623 -1.1200 0.0025
## sss 0.3806 1.0567 -1.6905 2.4517 0.7187
##
## ***********************
##
## Transformed coefficients:
## estimate standard.error lower.CI upper.CI p.value
## (Intercept) 15.3333 0.7472 13.8688 16.7978 0.0000
## amp 3.2138 1.0567 1.1427 5.2849 0.0024
## acr -0.1187 0.3288 -0.7631 0.5257 0.7181
# Draw the fitting curve 1h
## Draw the number of SE in a plot
circadianfigure <- ggplot.cosinor.lm(SEfit) + geom_bar(data = pSERGcircadian, aes(x = hourofSE, y = events), stat = 'identity', color = "#0076c0", fill = "#67771a", alpha = 0.8) +
scale_x_discrete(limits = c(1,24)) +
theme(panel.background = element_rect(fill = "white"),
axis.text=element_text(size = 18, color = "black", face = "bold"),
axis.title=element_text(size = 24, color = "black", face = "bold"),
axis.ticks.length = unit(0.4, "cm"), axis.ticks = element_line(size = 1.2)) +
scale_x_continuous(breaks = c(4, 8, 12, 16, 20, 24)) +
scale_y_continuous(breaks = c(5, 10, 15, 20, 25)) +
labs(x = "Time of day (hour)", y = "Number of rSE episodes") +
geom_segment(aes(x = 0, xend = 24, y = SEfit$coefficients[1], yend = SEfit$coefficients[1])) +
geom_segment(aes(x = 11.45, xend = 11.45, y = SEfit$coefficients[1], yend = (SEfit$coefficients[1] + SEfit$coefficients[2])), size=1.2) +
geom_segment(aes(x = 23.45, xend = 23.45, y = SEfit$coefficients[1], yend = (SEfit$coefficients[1] - SEfit$coefficients[2])), size=1.2)
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.
circadianfigure

# Interactive graph
ggplotly(circadianfigure)
#Draw the data and the fitting curve
circadianonly <- ggplot.cosinor.lm(SEfit) +
geom_point(data = circadiandata, aes(hourofSE, numberSE), color = "#67771a", size = 4) +
scale_x_discrete(limits = c(1,24)) +
theme(panel.background = element_rect(fill = "white"),
axis.text = element_text(size = 18, color = "black", face = "bold"),
axis.title = element_text(size = 24, color = "black", face = "bold"),
axis.ticks.length = unit(0.4, "cm"), axis.ticks = element_line(size = 1.2),
axis.line = element_line(color= "black", size = 1.2)) +
scale_x_continuous(breaks = c(4, 8, 12, 16, 20, 24)) +
labs(x = "Time of day (hour)", y = "Number of rSE episodes") +
geom_segment(aes(x = 0, xend = 24, y = SEfit$coefficients[1], yend = SEfit$coefficients[1])) +
geom_segment(aes(x = 10.55, xend = 10.55, y = SEfit$coefficients[1], yend = (SEfit$coefficients[1] + SEfit$coefficients[2]))) +
geom_segment(aes(x = 22.55, xend = 22.55, y = SEfit$coefficients[1], yend = (SEfit$coefficients[1] - SEfit$coefficients[2])))
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.
circadianonly

# Interactive graph
ggplotly(circadianonly)
Comparison patients with and without a neurologic history
## Proportion of patients with some neurological history per time bin
neurologichistory1h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 1, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 1, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 1, ]$none)[2])
neurologichistory2h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 2, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 2, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 2, ]$none)[2])
neurologichistory3h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 3, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 3, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 3, ]$none)[2])
neurologichistory4h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 4, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 4, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 4, ]$none)[2])
neurologichistory5h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 5, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 5, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 5, ]$none)[2])
neurologichistory6h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 6, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 6, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 6, ]$none)[2])
neurologichistory7h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 7, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 7, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 7, ]$none)[2])
neurologichistory8h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 8, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 8, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 8, ]$none)[2])
neurologichistory9h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 9, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 9, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 9, ]$none)[2])
neurologichistory10h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 10, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 10, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 10, ]$none)[2])
neurologichistory11h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 11, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 11, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 11, ]$none)[2])
neurologichistory12h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 12, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 12, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 12, ]$none)[2])
neurologichistory13h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 13, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 13, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 13, ]$none)[2])
neurologichistory14h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 14, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 14, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 14, ]$none)[2])
neurologichistory15h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 15, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 15, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 15, ]$none)[2])
neurologichistory16h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 16, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 16, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 16, ]$none)[2])
neurologichistory17h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 17, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 17, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 17, ]$none)[2])
neurologichistory18h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 18, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 18, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 18, ]$none)[2])
neurologichistory19h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 19, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 19, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 19, ]$none)[2])
neurologichistory20h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 20, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 20, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 20, ]$none)[2])
neurologichistory21h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 21, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 21, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 21, ]$none)[2])
neurologichistory22h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 22, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 22, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 22, ]$none)[2])
neurologichistory23h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 23, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 23, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 23, ]$none)[2])
neurologichistory24h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 24, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 24, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 24, ]$none)[2])
circadiandata$neurologic <- c(neurologichistory1h, neurologichistory2h, neurologichistory3h, neurologichistory4h, neurologichistory5h, neurologichistory6h, neurologichistory7h, neurologichistory8h, neurologichistory9h, neurologichistory10h, neurologichistory11h, neurologichistory12h, neurologichistory13h, neurologichistory14h, neurologichistory15h, neurologichistory16h, neurologichistory17h, neurologichistory18h, neurologichistory19h, neurologichistory20h, neurologichistory21h, neurologichistory22h, neurologichistory23h, neurologichistory24h)
circadiandata$neurologic <- ifelse(is.na(circadiandata$neurologic) == TRUE, 0, circadiandata$neurologic)
circadiandata$numberneurologic <- circadiandata$numberSE * circadiandata$neurologic
## Cosinor analysis
SEfitcomparisonneurologic <- cosinor.lm(numberSE ~ time(hourofSE) + numberneurologic + amp.acro(numberneurologic), data = circadiandata, period = 24)
summary(SEfitcomparisonneurologic)
## Raw model coefficients:
## estimate standard.error lower.CI upper.CI p.value
## (Intercept) 5.3261 1.6303 2.1307 8.5215 0.0011
## numberneurologic 0.9541 0.1597 0.6412 1.2671 0.0000
## rrr 0.0079 1.7118 -3.3471 3.3630 0.9963
## sss -1.4530 2.9021 -7.1411 4.2351 0.6166
## numberneurologic:rrr -0.2224 0.1635 -0.5427 0.0980 0.1737
## numberneurologic:sss 0.0512 0.2722 -0.4823 0.5846 0.8509
##
## ***********************
##
## Transformed coefficients:
## estimate standard.error lower.CI upper.CI
## (Intercept) 5.3261 1.6303 2.1307 8.5215
## [numberneurologic = 1] 0.9541 0.1597 0.6412 1.2671
## amp 1.4530 2.9012 -4.2332 7.1392
## [numberneurologic = 1]:amp 1.4181 2.6385 -3.7533 6.5896
## acr -1.5653 1.1793 -3.8766 0.7460
## [numberneurologic = 1]:acr 1.4190 1.0970 -0.7311 3.5691
## p.value
## (Intercept) 0.0011
## [numberneurologic = 1] 0.0000
## amp 0.6165
## [numberneurologic = 1]:amp 0.5909
## acr 0.1844
## [numberneurologic = 1]:acr 0.1958
test_cosinor(SEfitcomparisonneurologic, "numberneurologic", param = "amp")
## Global test:
## Statistic:
## [1] 0.01
##
##
## P-value:
## [1] 0.9234
##
## Individual tests:
## Statistic:
## [1] -0.1
##
##
## P-value:
## [1] 0.9234
##
## Estimate and confidence interval[1] "-0.03 (-0.75 to 0.68)"
test_cosinor(SEfitcomparisonneurologic, "numberneurologic", param = "acr")
## Global test:
## Statistic:
## [1] 86.37
##
##
## P-value:
## [1] 0
##
## Individual tests:
## Statistic:
## [1] 9.29
##
##
## P-value:
## [1] 0
##
## Estimate and confidence interval[1] "2.98 (2.35 to 3.61)"
##Graph
neurologichistoryfigure <- ggplot.cosinor.lm(SEfitcomparisonneurologic, x_str = "numberneurologic") +
geom_bar(data = pSERGcircadian[pSERGcircadian$none == 0, ], aes(x = hourofSE, y = events), stat = 'identity', color = "#5698a3", fill = "#5698a3", alpha = 0.2) +
geom_bar(data = pSERGcircadian[pSERGcircadian$none == 1, ], aes(x = hourofSE, y = events), stat = 'identity', color = "#a30234", fill = "#a30234", alpha = 0.2) +
scale_x_discrete(limits = c(1,24)) +
theme(panel.background=element_rect(fill = "white"),
axis.text=element_text(size = 18, color = "black", face = "bold"),
axis.title=element_text(size = 24, color = "black", face = "bold"),
axis.ticks.length=unit(0.4, "cm"), axis.ticks = element_line(size = 1.2),
legend.position = "none") +
scale_x_continuous(breaks = c(4, 8, 12, 16, 20, 24)) +
scale_y_continuous(breaks = c(5, 10, 15, 20)) +
labs(x = "Time of day (hour)", y = "Number of rSE episodes") +
geom_segment(aes(x = 0, xend = 24, y = SEfitcomparisonneurologic$coefficients[1], yend = SEfitcomparisonneurologic$coefficients[1]), color = "#a30234") +
geom_segment(aes(x = 18.2, xend = 18.2, y = SEfitcomparisonneurologic$coefficients[1], yend = (SEfitcomparisonneurologic$coefficients[1] + SEfitcomparisonneurologic$coefficients[3])), color = "#a30234", size=1.2) +
geom_segment(aes(x = 6.2, xend = 6.2, y = SEfitcomparisonneurologic$coefficients[1], yend = (SEfitcomparisonneurologic$coefficients[1] - SEfitcomparisonneurologic$coefficients[3])), color = "#a30234", size=1.2) +
geom_segment(aes(x = 0, xend = 24, y = SEfitcomparisonneurologic$coefficients[1] + SEfitcomparisonneurologic$coefficients[2], yend = SEfitcomparisonneurologic$coefficients[1] + SEfitcomparisonneurologic$coefficients[2]), color = "#5698a3") +
geom_segment(aes(x = 17.4, xend = 17.4, y = SEfitcomparisonneurologic$coefficients[1] + SEfitcomparisonneurologic$coefficients[2], yend = (SEfitcomparisonneurologic$coefficients[1] + SEfitcomparisonneurologic$coefficients[2] + SEfitcomparisonneurologic$coefficients[4])), color = "#5698a3", size=1.2) + geom_segment(aes(x = 5.4, xend = 5.4, y = SEfitcomparisonneurologic$coefficients[1] + SEfitcomparisonneurologic$coefficients[2], yend = (SEfitcomparisonneurologic$coefficients[1] + SEfitcomparisonneurologic$coefficients[2] - SEfitcomparisonneurologic$coefficients[4])), color = "#5698a3", size=1.2)
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.
neurologichistoryfigure

# Interactive graph
ggplotly(neurologichistoryfigure)
Comparison patients with and without a history of epilepsy
## Proportion of patients with a history of epilepsy per time bin
epilepsyhistory1h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 1, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 1, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 1, ]$priorepilepsy)[2])
epilepsyhistory2h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 2, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 2, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 2, ]$priorepilepsy)[2])
epilepsyhistory3h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 3, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 3, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 3, ]$priorepilepsy)[2])
epilepsyhistory4h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 4, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 4, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 4, ]$priorepilepsy)[2])
epilepsyhistory5h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 5, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 5, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 5, ]$priorepilepsy)[2])
epilepsyhistory6h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 6, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 6, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 6, ]$priorepilepsy)[2])
epilepsyhistory7h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 7, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 7, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 7, ]$priorepilepsy)[2])
epilepsyhistory8h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 8, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 8, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 8, ]$priorepilepsy)[2])
epilepsyhistory9h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 9, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 9, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 9, ]$priorepilepsy)[2])
epilepsyhistory10h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 10, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 10, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 10, ]$priorepilepsy)[2])
epilepsyhistory11h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 11, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 11, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 11, ]$priorepilepsy)[2])
epilepsyhistory12h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 12, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 12, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 12, ]$priorepilepsy)[2])
epilepsyhistory13h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 13, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 13, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 13, ]$priorepilepsy)[2])
epilepsyhistory14h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 14, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 14, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 14, ]$priorepilepsy)[2])
epilepsyhistory15h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 15, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 15, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 15, ]$priorepilepsy)[2])
epilepsyhistory16h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 16, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 16, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 16, ]$priorepilepsy)[2])
epilepsyhistory17h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 17, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 17, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 17, ]$priorepilepsy)[2])
epilepsyhistory18h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 18, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 18, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 18, ]$priorepilepsy)[2])
epilepsyhistory19h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 19, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 19, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 19, ]$priorepilepsy)[2])
epilepsyhistory20h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 20, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 20, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 20, ]$priorepilepsy)[2])
epilepsyhistory21h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 21, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 21, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 21, ]$priorepilepsy)[2])
epilepsyhistory22h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 22, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 22, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 22, ]$priorepilepsy)[2])
epilepsyhistory23h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 23, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 23, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 23, ]$priorepilepsy)[2])
epilepsyhistory24h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 24, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 24, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 24, ]$priorepilepsy)[2])
circadiandata$epilepsy <- c(epilepsyhistory1h, epilepsyhistory2h, epilepsyhistory3h, epilepsyhistory4h, epilepsyhistory5h, epilepsyhistory6h, epilepsyhistory7h, epilepsyhistory8h, epilepsyhistory9h, epilepsyhistory10h, epilepsyhistory11h, epilepsyhistory12h, epilepsyhistory13h, epilepsyhistory14h, epilepsyhistory15h, epilepsyhistory16h, epilepsyhistory17h, epilepsyhistory18h, epilepsyhistory19h, epilepsyhistory20h, epilepsyhistory21h, epilepsyhistory22h, epilepsyhistory23h, epilepsyhistory24h)
circadiandata$epilepsy <- ifelse(is.na(circadiandata$epilepsy) == TRUE, 0, circadiandata$epilepsy)
circadiandata$numberepilepsy <- circadiandata$numberSE * circadiandata$epilepsy
## Cosinor analysis
SEfitcomparisonepilepsy <- cosinor.lm(numberSE ~ time(hourofSE) + numberepilepsy + amp.acro(numberepilepsy), data = circadiandata, period = 24)
summary(SEfitcomparisonepilepsy)
## Raw model coefficients:
## estimate standard.error lower.CI upper.CI p.value
## (Intercept) 7.4248 1.4394 4.6037 10.2460 0.0000
## numberepilepsy 1.0437 0.2000 0.6518 1.4357 0.0000
## rrr -1.3100 1.4850 -4.2206 1.6006 0.3777
## sss -0.6892 2.7395 -6.0585 4.6800 0.8014
## numberepilepsy:rrr -0.1829 0.1869 -0.5492 0.1834 0.3277
## numberepilepsy:sss -0.0757 0.3557 -0.7728 0.6213 0.8314
##
## ***********************
##
## Transformed coefficients:
## estimate standard.error lower.CI upper.CI p.value
## (Intercept) 7.4248 1.4394 4.6037 10.2460 0.0000
## [numberepilepsy = 1] 1.0437 0.2000 0.6518 1.4357 0.0000
## amp 1.4803 1.7568 -1.9630 4.9235 0.3994
## [numberepilepsy = 1]:amp 1.6775 1.5324 -1.3260 4.6810 0.2737
## acr 0.4843 1.7386 -2.9233 3.8920 0.7806
## [numberepilepsy = 1]:acr 0.4735 1.3506 -2.1735 3.1206 0.7259
test_cosinor(SEfitcomparisonepilepsy, "numberepilepsy", param = "amp")
## Global test:
## Statistic:
## [1] 0.66
##
##
## P-value:
## [1] 0.4182
##
## Individual tests:
## Statistic:
## [1] 0.81
##
##
## P-value:
## [1] 0.4182
##
## Estimate and confidence interval[1] "0.2 (-0.28 to 0.67)"
test_cosinor(SEfitcomparisonepilepsy, "numberepilepsy", param = "acr")
## Global test:
## Statistic:
## [1] 0
##
##
## P-value:
## [1] 0.978
##
## Individual tests:
## Statistic:
## [1] -0.03
##
##
## P-value:
## [1] 0.978
##
## Estimate and confidence interval[1] "-0.01 (-0.78 to 0.76)"
##Graph
epilepsyhistoryfigure <- ggplot.cosinor.lm(SEfitcomparisonepilepsy, x_str = "numberepilepsy") +
geom_bar(data = pSERGcircadian[pSERGcircadian$priorepilepsy == 1, ], aes(x = hourofSE, y = events), stat = 'identity', color = "#5698a3", fill = "#5698a3", alpha = 0.2) +
geom_bar(data = pSERGcircadian[pSERGcircadian$priorepilepsy == 0, ], aes(x = hourofSE, y = events), stat = 'identity', color = "#a30234", fill = "#a30234", alpha = 0.2) +
scale_x_discrete(limits = c(1,24)) +
theme(panel.background=element_rect(fill = "white"),
axis.text=element_text(size = 18, color = "black", face = "bold"),
axis.title=element_text(size = 24, color = "black", face = "bold"),
axis.ticks.length=unit(0.4, "cm"), axis.ticks = element_line(size = 1.2),
legend.position = "none") +
scale_x_continuous(breaks = c(4, 8, 12, 16, 20, 24)) +
scale_y_continuous(breaks = c(5, 10, 15)) +
labs(x = "Time of day (hour)", y = "Number of rSE episodes") +
geom_segment(aes(x = 0, xend = 24, y = SEfitcomparisonepilepsy$coefficients[1], yend = SEfitcomparisonepilepsy$coefficients[1]), color = "#a30234") +
geom_segment(aes(x = 13.9, xend = 13.9, y = SEfitcomparisonepilepsy$coefficients[1], yend = (SEfitcomparisonepilepsy$coefficients[1] + SEfitcomparisonepilepsy$coefficients[3])), color = "#a30234", size=1.2) +
geom_segment(aes(x = 1.9, xend = 1.9, y = SEfitcomparisonepilepsy$coefficients[1], yend = (SEfitcomparisonepilepsy$coefficients[1] - SEfitcomparisonepilepsy$coefficients[3])), color = "#a30234", size=1.2) +
geom_segment(aes(x = 0, xend = 24, y = SEfitcomparisonepilepsy$coefficients[1] + SEfitcomparisonepilepsy$coefficients[2], yend = SEfitcomparisonepilepsy$coefficients[1] + SEfitcomparisonepilepsy$coefficients[2]), color = "#5698a3") +
geom_segment(aes(x = 13.8, xend = 13.8, y = SEfitcomparisonepilepsy$coefficients[1] + SEfitcomparisonepilepsy$coefficients[2], yend = (SEfitcomparisonepilepsy$coefficients[1] + SEfitcomparisonepilepsy$coefficients[2] + SEfitcomparisonepilepsy$coefficients[4])), color = "#5698a3", size=1.2) + geom_segment(aes(x = 1.8, xend = 1.8, y = SEfitcomparisonepilepsy$coefficients[1] + SEfitcomparisonepilepsy$coefficients[2], yend = (SEfitcomparisonepilepsy$coefficients[1] + SEfitcomparisonepilepsy$coefficients[2] - SEfitcomparisonepilepsy$coefficients[4])), color = "#5698a3", size=1.2)
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.
epilepsyhistoryfigure

# Interactive graph
ggplotly(epilepsyhistoryfigure)
Comparison patients with SE onset in and out of the hospital
## Proportion of patients with onset of SE out of the hospital
outofhospitalonset1h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 1, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 1, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 1, ]$HOSPITALONSET)[2])
outofhospitalonset2h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 2, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 2, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 2, ]$HOSPITALONSET)[2])
outofhospitalonset3h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 3, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 3, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 3, ]$HOSPITALONSET)[2])
outofhospitalonset4h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 4, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 4, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 4, ]$HOSPITALONSET)[2])
outofhospitalonset5h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 5, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 5, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 5, ]$HOSPITALONSET)[2])
outofhospitalonset6h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 6, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 6, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 6, ]$HOSPITALONSET)[2])
outofhospitalonset7h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 7, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 7, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 7, ]$HOSPITALONSET)[2])
outofhospitalonset8h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 8, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 8, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 8, ]$HOSPITALONSET)[2])
outofhospitalonset9h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 9, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 9, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 9, ]$HOSPITALONSET)[2])
outofhospitalonset10h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 10, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 10, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 10, ]$HOSPITALONSET)[2])
outofhospitalonset11h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 11, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 11, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 11, ]$HOSPITALONSET)[2])
outofhospitalonset12h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 12, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 12, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 12, ]$HOSPITALONSET)[2])
outofhospitalonset13h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 13, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 13, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 13, ]$HOSPITALONSET)[2])
outofhospitalonset14h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 14, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 14, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 14, ]$HOSPITALONSET)[2])
outofhospitalonset15h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 15, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 15, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 15, ]$HOSPITALONSET)[2])
outofhospitalonset16h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 16, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 16, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 16, ]$HOSPITALONSET)[2])
outofhospitalonset17h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 17, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 17, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 17, ]$HOSPITALONSET)[2])
outofhospitalonset18h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 18, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 18, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 18, ]$HOSPITALONSET)[2])
outofhospitalonset19h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 19, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 19, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 19, ]$HOSPITALONSET)[2])
outofhospitalonset20h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 20, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 20, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 20, ]$HOSPITALONSET)[2])
outofhospitalonset21h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 21, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 21, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 21, ]$HOSPITALONSET)[2])
outofhospitalonset22h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 22, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 22, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 22, ]$HOSPITALONSET)[2])
outofhospitalonset23h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 23, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 23, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 23, ]$HOSPITALONSET)[2])
outofhospitalonset24h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 24, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 24, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 24, ]$HOSPITALONSET)[2])
circadiandata$outofhospitalonset <- c(outofhospitalonset1h, outofhospitalonset2h, outofhospitalonset3h,
outofhospitalonset4h, outofhospitalonset5h, outofhospitalonset6h,
outofhospitalonset7h, outofhospitalonset8h, outofhospitalonset9h,
outofhospitalonset10h, outofhospitalonset11h, outofhospitalonset12h,
outofhospitalonset13h, outofhospitalonset14h, outofhospitalonset15h,
outofhospitalonset16h, outofhospitalonset17h, outofhospitalonset18h,
outofhospitalonset19h, outofhospitalonset20h, outofhospitalonset21h,
outofhospitalonset22h, outofhospitalonset23h, outofhospitalonset24h)
circadiandata$outofhospitalonset <- ifelse(is.na(circadiandata$outofhospitalonset) == TRUE, 0, circadiandata$outofhospitalonset)
circadiandata$numberoutofhospitalonset <- circadiandata$numberSE * circadiandata$outofhospitalonset
## Cosinor analysis
SEfitcomparisonhospitalonset <- cosinor.lm(numberSE ~ time(hourofSE) + numberoutofhospitalonset + amp.acro(numberoutofhospitalonset), data = circadiandata, period = 24)
summary(SEfitcomparisonhospitalonset)
## Raw model coefficients:
## estimate standard.error lower.CI upper.CI
## (Intercept) 6.9350 3.4894 0.0958 13.7741
## numberoutofhospitalonset 0.7835 0.3002 0.1950 1.3719
## rrr -1.6279 4.2202 -9.8993 6.6436
## sss -2.3001 4.4855 -11.0916 6.4914
## numberoutofhospitalonset:rrr 0.0168 0.3766 -0.7213 0.7549
## numberoutofhospitalonset:sss 0.2380 0.3821 -0.5109 0.9869
## p.value
## (Intercept) 0.0469
## numberoutofhospitalonset 0.0091
## rrr 0.6997
## sss 0.6081
## numberoutofhospitalonset:rrr 0.9645
## numberoutofhospitalonset:sss 0.5334
##
## ***********************
##
## Transformed coefficients:
## estimate standard.error lower.CI
## (Intercept) 6.9350 3.4894 0.0958
## [numberoutofhospitalonset = 1] 0.7835 0.3002 0.1950
## amp 2.8179 4.1158 -5.2490
## [numberoutofhospitalonset = 1]:amp 2.6169 3.7532 -4.7393
## acr 0.9549 1.6259 -2.2317
## [numberoutofhospitalonset = 1]:acr 0.9076 1.6069 -2.2419
## upper.CI p.value
## (Intercept) 13.7741 0.0469
## [numberoutofhospitalonset = 1] 1.3719 0.0091
## amp 10.8848 0.4936
## [numberoutofhospitalonset = 1]:amp 9.9730 0.4857
## acr 4.1415 0.5570
## [numberoutofhospitalonset = 1]:acr 4.0570 0.5722
test_cosinor(SEfitcomparisonhospitalonset, "numberoutofhospitalonset", param = "amp")
## Global test:
## Statistic:
## [1] 0.23
##
##
## P-value:
## [1] 0.6307
##
## Individual tests:
## Statistic:
## [1] -0.48
##
##
## P-value:
## [1] 0.6307
##
## Estimate and confidence interval[1] "-0.2 (-1.02 to 0.62)"
test_cosinor(SEfitcomparisonhospitalonset, "numberoutofhospitalonset", param = "acr")
## Global test:
## Statistic:
## [1] 0.39
##
##
## P-value:
## [1] 0.5311
##
## Individual tests:
## Statistic:
## [1] -0.63
##
##
## P-value:
## [1] 0.5311
##
## Estimate and confidence interval[1] "-0.05 (-0.2 to 0.1)"
##Graph
hospitalonsetfigure <- ggplot.cosinor.lm(SEfitcomparisonhospitalonset, x_str = "numberoutofhospitalonset") +
geom_bar(data = pSERGcircadian[pSERGcircadian$HOSPITALONSET == "no", ], aes(x = hourofSE, y = events), stat = 'identity', color = "#5698a3", fill = "#5698a3", alpha = 0.2) +
geom_bar(data = pSERGcircadian[pSERGcircadian$HOSPITALONSET == "yes", ], aes(x = hourofSE, y = events), stat = 'identity', color = "#a30234", fill = "#a30234", alpha = 0.2) +
scale_x_discrete(limits = c(1,24)) +
theme(panel.background=element_rect(fill = "white"),
axis.text=element_text(size = 18, color = "black", face = "bold"),
axis.title=element_text(size = 24, color = "black", face = "bold"),
axis.ticks.length=unit(0.4, "cm"), axis.ticks = element_line(size = 1.2),
legend.position = "none") +
scale_x_continuous(breaks = c(4, 8, 12, 16, 20, 24)) +
scale_y_continuous(breaks = c(5, 10, 15)) +
labs(x = "Time of day (hour)", y = "Number of rSE episodes") +
geom_segment(aes(x = 0, xend = 24, y = SEfitcomparisonhospitalonset$coefficients[1], yend = SEfitcomparisonhospitalonset$coefficients[1]), color = "#a30234") +
geom_segment(aes(x = 15.7, xend = 15.7, y = SEfitcomparisonhospitalonset$coefficients[1], yend = (SEfitcomparisonhospitalonset$coefficients[1] + SEfitcomparisonhospitalonset$coefficients[3])), color = "#a30234", size=1.2) +
geom_segment(aes(x = 3.7, xend = 3.7, y = SEfitcomparisonhospitalonset$coefficients[1], yend = (SEfitcomparisonhospitalonset$coefficients[1] - SEfitcomparisonhospitalonset$coefficients[3])), color = "#a30234", size=1.2) +
geom_segment(aes(x = 0, xend = 24, y = SEfitcomparisonhospitalonset$coefficients[1] + SEfitcomparisonhospitalonset$coefficients[2], yend = SEfitcomparisonhospitalonset$coefficients[1] + SEfitcomparisonhospitalonset$coefficients[2]), color = "#5698a3") +
geom_segment(aes(x = 15.4, xend = 15.4, y = SEfitcomparisonhospitalonset$coefficients[1] + SEfitcomparisonhospitalonset$coefficients[2], yend = (SEfitcomparisonhospitalonset$coefficients[1] + SEfitcomparisonhospitalonset$coefficients[2] + SEfitcomparisonhospitalonset$coefficients[4])), color = "#5698a3", size=1.2) + geom_segment(aes(x = 3.4, xend = 3.4, y = SEfitcomparisonhospitalonset$coefficients[1] + SEfitcomparisonhospitalonset$coefficients[2], yend = (SEfitcomparisonhospitalonset$coefficients[1] + SEfitcomparisonhospitalonset$coefficients[2] - SEfitcomparisonhospitalonset$coefficients[4])), color = "#5698a3", size=1.2)
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.
hospitalonsetfigure

# Interactive graph
ggplotly(hospitalonsetfigure)
Time to first BZD
##Transform hours of SE to factor
pSERGcircadian$hourofSE<- as.factor(pSERGcircadian$hourofSE)
# Transform time to first BZD time to numeric
pSERGcircadian$BZDTIME.0 <- as.numeric(as.character(pSERGcircadian$BZDTIME.0))
## Warning: NAs introduced by coercion
#Summary data
summary(pSERGcircadian$BZDTIME.0)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.00 5.00 16.50 70.06 45.00 1605.00 28
######################Transform hours of SE to numeric############
pSERGcircadian$hourofSE <- as.numeric(pSERGcircadian$hourofSE)
#Cosinor analysis
BZDfit <- cosinor.lm(BZDTIME.0 ~ time(hourofSE), data = pSERGcircadian[!is.na(pSERGcircadian$hourofSE), ], period = 24)
summary(BZDfit)
## Raw model coefficients:
## estimate standard.error lower.CI upper.CI p.value
## (Intercept) 71.1615 10.2861 51.0012 91.3218 0.0000
## rrr 10.3922 14.6235 -18.2693 39.0536 0.4773
## sss 0.4276 14.3476 -27.6932 28.5483 0.9762
##
## ***********************
##
## Transformed coefficients:
## estimate standard.error lower.CI upper.CI p.value
## (Intercept) 71.1615 10.2861 51.0012 91.3218 0.0000
## amp 10.4010 14.6555 -18.3233 39.1252 0.4779
## acr 0.0411 1.3763 -2.6564 2.7386 0.9762
##Draw the number of SE in a plot
################Transform hours of SE to factor######
pSERGcircadian$hourofSE <- as.factor(pSERGcircadian$hourofSE)
timeBZD <- ggplot.cosinor.lm(BZDfit) +
geom_boxplot(data = pSERGcircadian[!is.na(pSERGcircadian$hourofSE), ], aes(x = hourofSE, y = BZDTIME.0), color = "#0076c0", fill = "#67771a", alpha = 0.5) +
scale_x_discrete(limits = c(1,24)) +
theme(panel.background=element_rect(fill = "white"),
axis.text=element_text(size = 16, color = "black", face = "bold"),
axis.title=element_text(size = 18, color = "black", face = "bold"),
axis.ticks.length = unit(0.4, "cm"), axis.ticks = element_line(size = 1.2)) +
scale_x_discrete(breaks = c("4", "8", "12", "16", "20", "24")) +
scale_y_continuous(breaks = c(0, 50, 100, 150, 200, 250, 300, 350)) +
coord_cartesian(ylim = c(0, 200)) +
labs(x = "Time of day (hour)", y= "Time to the first BZD (min)") +
geom_segment(aes(x = 0,xend = 24, y= BZDfit$coefficients[1], yend = BZDfit$coefficients[1])) +
geom_segment(aes(x = 0.5, xend = 0.5, y = BZDfit$coefficients[1], yend = BZDfit$coefficients[1] + BZDfit$coefficients[2]), size=1.2) +
geom_segment(aes(x = 12.5, xend = 12.5, y = BZDfit$coefficients[1], yend = BZDfit$coefficients[1] - BZDfit$coefficients[2]), size=1.2)
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.
timeBZD
## Warning: Removed 28 rows containing non-finite values (stat_boxplot).

# Interactive graph
ggplotly(timeBZD)
## Warning: Removed 28 rows containing non-finite values (stat_boxplot).
Time to first non-BZD AED
## Transform hours of SE to factor
pSERGcircadian$hourofSE <- as.factor(pSERGcircadian$hourofSE)
# Transform time to first AED time to numeric
pSERGcircadian$AEDTIME.0 <- as.numeric(as.character(pSERGcircadian$AEDTIME.0))
## Warning: NAs introduced by coercion
# Summary data
summary(pSERGcircadian$AEDTIME.0)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.0 35.0 66.0 165.4 150.0 4320.0 27
## Transform hours of SE to numeric
pSERGcircadian$hourofSE <- as.numeric(pSERGcircadian$hourofSE)
# Cosinor analysis
ASMfit <- cosinor.lm(AEDTIME.0~time(hourofSE), data=pSERGcircadian, period=24)
summary(ASMfit)
## Raw model coefficients:
## estimate standard.error lower.CI upper.CI p.value
## (Intercept) 171.9320 19.5132 133.6869 210.1771 0.0000
## rrr 50.7544 28.0200 -4.1639 105.6727 0.0701
## sss 44.6559 26.9189 -8.1042 97.4159 0.0971
##
## ***********************
##
## Transformed coefficients:
## estimate standard.error lower.CI upper.CI p.value
## (Intercept) 171.9320 19.5132 133.6869 210.1771 0.0000
## amp 67.6029 28.5555 11.6352 123.5706 0.0179
## acr 0.7216 0.3898 -0.0424 1.4855 0.0641
################Transform hours of SE to factor######
pSERGcircadian$hourofSE <- as.factor(pSERGcircadian$hourofSE)
#Draw the data and the fitting curve
timeASM <- ggplot.cosinor.lm(ASMfit) +
geom_boxplot(data = pSERGcircadian[!is.na(pSERGcircadian$hourofSE), ], aes(hourofSE, AEDTIME.0), color = "#0076c0", fill = "#67771a", alpha = 0.5) +
scale_x_discrete(limits = c(1,24)) +
scale_y_discrete(limits = c(0,800)) +
theme(panel.background = element_rect(fill = "white"),
axis.text = element_text(size = 16, color = "black", face = "bold"),
axis.title = element_text(size = 18, color = "black", face = "bold"),
axis.ticks.length = unit(0.4, "cm"), axis.ticks = element_line(size = 1.2)) +
scale_x_discrete(breaks = c("4", "8", "12", "16", "20", "24")) +
scale_y_continuous(breaks = c(0, 100, 200, 300, 400, 500, 600)) +
coord_cartesian(ylim = c(0, 600)) +
labs(x = "Time of day (hour)", y = "Time to the first non-BZD AED (min)") +
geom_segment(aes(x = 0, xend = 24, y = ASMfit$coefficients[1], yend = ASMfit$coefficients[1])) +
geom_segment(aes(x = 2.75, xend = 2.75, y = ASMfit$coefficients[1], yend = ASMfit$coefficients[1] + ASMfit$coefficients[2]), size=1.2) +
geom_segment(aes(x = 14.75, xend = 14.75, y = ASMfit$coefficients[1], yend = ASMfit$coefficients[1] - ASMfit$coefficients[2]), size=1.2)
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.
timeASM
## Warning: Removed 27 rows containing non-finite values (stat_boxplot).

# Interactive graph
ggplotly(timeASM)
## Warning: Removed 27 rows containing non-finite values (stat_boxplot).
Alternative hour definition
#################Create variable sensitivity analysis hour of SE######################
## Define the newly defined hour
pSERGcircadian$minutes <- pSERGcircadian$TIMESEIZURE_MINUTES
pSERGcircadian$minutesashours <- pSERGcircadian$minutes / 60
pSERGcircadian$minutesashours <- ifelse(is.na(pSERGcircadian$minutesashours) == TRUE | pSERGcircadian$minutesashours > 1, 0, pSERGcircadian$minutesashours)
## Create variable new hour
pSERGcircadian$newhour <- as.numeric(pSERGcircadian$hourofSE) + pSERGcircadian$minutesashours
## Rename hours
pSERGcircadian$newh[pSERGcircadian$newhour >= 24.5 | pSERGcircadian$newhour < 1.5] <- 1
pSERGcircadian$newh[pSERGcircadian$newhour >= 1.5 & pSERGcircadian$newhour < 2.5] <- 2
pSERGcircadian$newh[pSERGcircadian$newhour >= 2.5 & pSERGcircadian$newhour < 3.5] <- 3
pSERGcircadian$newh[pSERGcircadian$newhour >= 3.5 & pSERGcircadian$newhour < 4.5] <- 4
pSERGcircadian$newh[pSERGcircadian$newhour >= 4.5 & pSERGcircadian$newhour < 5.5] <- 5
pSERGcircadian$newh[pSERGcircadian$newhour >= 5.5 & pSERGcircadian$newhour < 6.5] <- 6
pSERGcircadian$newh[pSERGcircadian$newhour >= 6.5 & pSERGcircadian$newhour < 7.5] <- 7
pSERGcircadian$newh[pSERGcircadian$newhour >= 7.5 & pSERGcircadian$newhour < 8.5] <- 8
pSERGcircadian$newh[pSERGcircadian$newhour >= 8.5 & pSERGcircadian$newhour < 9.5] <- 9
pSERGcircadian$newh[pSERGcircadian$newhour >= 9.5 & pSERGcircadian$newhour < 10.5] <- 10
pSERGcircadian$newh[pSERGcircadian$newhour >= 10.5 & pSERGcircadian$newhour < 11.5] <- 11
pSERGcircadian$newh[pSERGcircadian$newhour >= 11.5 & pSERGcircadian$newhour < 12.5] <- 12
pSERGcircadian$newh[pSERGcircadian$newhour >= 12.5 & pSERGcircadian$newhour < 13.5] <- 13
pSERGcircadian$newh[pSERGcircadian$newhour >= 13.5 & pSERGcircadian$newhour < 14.5] <- 14
pSERGcircadian$newh[pSERGcircadian$newhour >= 14.5 & pSERGcircadian$newhour < 15.5] <- 15
pSERGcircadian$newh[pSERGcircadian$newhour >= 15.5 & pSERGcircadian$newhour < 16.5] <- 16
pSERGcircadian$newh[pSERGcircadian$newhour >= 16.5 & pSERGcircadian$newhour < 17.5] <- 17
pSERGcircadian$newh[pSERGcircadian$newhour >= 17.5 & pSERGcircadian$newhour < 18.5] <- 18
pSERGcircadian$newh[pSERGcircadian$newhour >= 18.5 & pSERGcircadian$newhour < 19.5] <- 19
pSERGcircadian$newh[pSERGcircadian$newhour >= 19.5 & pSERGcircadian$newhour < 20.5] <- 20
pSERGcircadian$newh[pSERGcircadian$newhour >= 20.5 & pSERGcircadian$newhour < 21.5] <- 21
pSERGcircadian$newh[pSERGcircadian$newhour >= 21.5 & pSERGcircadian$newhour < 22.5] <- 22
pSERGcircadian$newh[pSERGcircadian$newhour >= 22.5 & pSERGcircadian$newhour < 23.5] <- 23
pSERGcircadian$newh[pSERGcircadian$newhour >= 23.5 & pSERGcircadian$newhour < 24.5] <- 24
# Test (Kolmogorov-Smirnov) if the distribution of SE during the day is consistent with a uniform distribution
ks.test(pSERGcircadian$newh, "punif", 0, 24)
## Warning in ks.test(pSERGcircadian$newh, "punif", 0, 24): ties should not be
## present for the Kolmogorov-Smirnov test
##
## One-sample Kolmogorov-Smirnov test
##
## data: pSERGcircadian$newh
## D = 0.082428, p-value = 0.01347
## alternative hypothesis: two-sided
#####################################################################################
# Create another database with the circadian data
newoccurrences <- count(pSERGcircadian, "newh")
# Make it a dataframe
newcircadiandata <- data.frame(newoccurrences)
newcircadiandata$numberSE <- newcircadiandata$freq
newcircadiandata <- newcircadiandata[ , c("newh", "numberSE")]
#####################################################################################
# Cosinor analysis
newSEfit <- cosinor.lm(numberSE ~ time(newh), data = newcircadiandata, period = 24)
summary(newSEfit)
## Raw model coefficients:
## estimate standard.error lower.CI upper.CI p.value
## (Intercept) 15.3333 0.7375 13.8878 16.7788 0.0000
## rrr -3.4051 1.0430 -5.4493 -1.3608 0.0011
## sss -0.0730 1.0430 -2.1172 1.9712 0.9442
##
## ***********************
##
## Transformed coefficients:
## estimate standard.error lower.CI upper.CI p.value
## (Intercept) 15.3333 0.7375 13.8878 16.7788 0.0000
## amp 3.4058 1.0430 1.3616 5.4501 0.0011
## acr 0.0214 0.3062 -0.5788 0.6216 0.9442
################Every hour new definition ###############
## Draw the number of SE episodes in a plot with stacked units
newdistributionhour <- ggplot(pSERGcircadian, aes(x = round(newh), y = events)) + geom_bar(stat = 'identity', color = "#0076c0", fill = "#67771a") +
scale_x_discrete(limits = c(1,24)) +
coord_cartesian(xlim = c(1, 24)) +
theme(panel.background = element_rect(fill = "white"),
axis.text = element_text(size = 18, color = "black", face = "bold"),
axis.title = element_text(size = 24, color = "black", face = "bold"),
axis.ticks.length = unit(0.4, "cm"), axis.ticks = element_line(size = 1.2)) +
scale_x_continuous(breaks = c(4, 8, 12, 16, 20, 24)) +
scale_y_continuous(breaks = c(5, 10, 15, 20, 25)) +
labs(x= "Time of day (hour)", y= "Number of rSE episodes")
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.
newdistributionhour

# Interactive graph
ggplotly(newdistributionhour)
# Draw the fitting curve 1h
## Draw the number of SE in a plot with the redefined hour
newcircadianfigure <- ggplot.cosinor.lm(SEfit) + geom_bar(data = pSERGcircadian, aes(x = newh, y = events), stat = 'identity', color = "#0076c0", fill = "#67771a", alpha = 0.8) +
scale_x_discrete(limits = c(1,24)) +
theme(panel.background = element_rect(fill = "white"),
axis.text=element_text(size = 18, color = "black", face = "bold"),
axis.title=element_text(size = 24, color = "black", face = "bold"),
axis.ticks.length = unit(0.4, "cm"), axis.ticks = element_line(size = 1.2)) +
scale_x_continuous(breaks = c(4, 8, 12, 16, 20, 24)) +
scale_y_continuous(breaks = c(5, 10, 15, 20, 25)) +
labs(x = "Time of day (hour)", y = "Number of rSE episodes") +
geom_segment(aes(x = 0, xend = 24, y = SEfit$coefficients[1], yend = SEfit$coefficients[1])) +
geom_segment(aes(x = 11.7, xend = 11.7, y = SEfit$coefficients[1], yend = (SEfit$coefficients[1] + SEfit$coefficients[2])), size=1.2) +
geom_segment(aes(x = 23.7, xend = 23.7, y = SEfit$coefficients[1], yend = (SEfit$coefficients[1] - SEfit$coefficients[2])), size=1.2)
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.
newcircadianfigure

# Interactive graph
ggplotly(newcircadianfigure)